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tions between cells in the form of a Fokker-Planck equation for the evolution of the cell probability 
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I. INTRODUCTION 

Biological cell dynamics has been studied at two main scales of description. The macroscopic level provides one 
with a coarse-grained treatment of biological cells through their macroscopically averaged quantities such as local 
density of cells [H E IS 

The macroscopic scale is large in comparison with the typical size of a cell. Macroscopic 
models are usually continuous and utilize families of differential or integro-differential equations to describe "fields" 
of interaction. A much more detailed approach is needed at the second, microscopic level which takes into account 
stochastic fluctuations of the shape of each individual cell. 

Discrete models describe individual (microscopic) behaviors of cells. They are often applied to microscale events 
where a small number of elements can have a large (and stochastic) impact on a system. For example, while many 
periodic growth patterns can be modeled using continuous methods, patterns which depend sensitively on interaction 
between cells and substrate are best modeled with discrete methods. Simplest discrete models describe cells as point- 
wise objects. Some bacteria are self-propelled and do not change considerably their shape during motion (e.g. E. 
Coli [SlU and M. xanthus JEHbacteria) . They can be successfully represented as point-wise objects undergoing 
reorientation while moving In contrast, some other bacteria (e.g. Dictyostelium discoideum 0]) experience 

essential random fluctuations of their shapes and need to be treated as extended objects of variable shapes. 

One of the microscopic models dealing with differential adhesion and shape fluctuations is a Cellular Potts Model 
(CPM) which is an extension of the well known Potts Model from statistical mechanics 0,0]. In this model each 
biological cell is represented by a cluster of pixels (spins). The CPM has been used to simulate various biolo gica l 
phenomena such as cell sorting |llL , fruiting body formation of the slime mold Dictyostelium discoideum [Till 14) , 
angiogenesis |l5j . cancer invasion |l6j. chondrogenesis in embryonic vertebrate limbs |l7l Il8| . and many others. 
(Different applications of the CPM have been reviewed in 0.) Recently a new alternative model was suggested [2(j 
which represents a cell as collection of subcellular elements which interact with each other through phenomenological 
; i intra- and intercellular potentials. 

In addition to short range cell-cell adhesion and interactions between cells and their surrounding extracellular matrix 
(haptotaxis), cell interact at long range through signal transmission and reception mediated by a diffusing chemical 
field (chemotaxis). Continuous macroscopic Keller-Segel model of the evolution of the density of cells with chemotactic 
interactions have been extensively studied Q, 0, S U over the years. In particular, it has been successfully applied 
to the description of Escherichia coli bacteria aggregation due to chemotaxis in |2|. The drawback of continuous 
models is that they have a lower resolution than discrete models. However, their advantage is the availability of a 
large set of analytical and numerical tools for analyzing solutions of the corresponding nonlinear partial differential 
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equations (PDEs). By contrast, the analytical study of discrete models is often impossibly complicated, and their 
computational implementation is often much less efficient in comparison with numerical methods available for PDEs. 
It is thus important, for numerical, analytical, as well as conceptual reasons to establish connections between various 
discrete and continuous models of the same biological problem. 

There is a vast literature on studying continuous limits of point-wise discrete microscopic models. In particular, 
classical Keller-Segel model has been derived from a model with point-wise representation for cells undergoing random 
walk |sl l2lL |22^| . However, much less work has been done on deriving macroscopic limits of microscopic models which 
treat cells as extended objects. One of the first attempts at combining microscopic and macroscopic levels of description 
of cellular dynamics has been described in [2j| where the diffusion coefficient for a collection of noninteracting randomly 
moving cells has been derived from a one dimensional CPM. Recently a microscopic limit of subcellular elements model 
|20| was derived in the form of continuous advection-diffusion equation for cellular density. In the present paper, we 
establish a connection between a one-dimensional CPM of a cell moving on a substrate and reacting to a chemical 
field, and a Fokker-Planck equation for the cell probability density function. This equation is then reduced to the 
classical macroscopic Keller-Segel equation. In particular, we derive all coefficients of the Keller-Segel model from 
parameters of the CPM. We also compare Monte Carlo simulations for the CPM with numerics for the Keller-Segel 
model to support our theoretical results. 

Unified multiscale approach, described in this paper and based on combining microscopic and macroscopic models, 
can be applied to studying such biological phenomena as streaming in Dictyostelium discoideum. In starved pop- 
ulations of Dictyostelium amoebae, cells produce and detect a communication chemical (cAMP). The movement of 
Dictyostelium cells changes from a random walk to a directed walk up the cAMP gradient resulting in formation of 
streams of cells towards the aggregation center (see Fig. and subsequent formation of multi-cellular fruiting body. 
Figure ^3 shows cells' movement from left to right in response to waves of cAMP travelling through the aggregation 
stream from right to left. The cAMP gradient on the up-down direction is very small and could be ig nored. Figure^ 
schematically demonstrates the main features of the cell movement. Unlike differential adhesion |lllll2]. chemotactic 
cell motion is highly organized over a length scale significantly larger than the size of a single cell. (For details about 
modeling Dictyostelium discoideum fruiting body formation see e.g. p^.ll4ll2^.l25|). 

The paper is organized as follows. In Section [H] we describe a one dimensional CPM with chemotaxis. In Section 
IIIII we derive from the Monte Carlo dynamics of the CPM the discrete master equation for the probability density 
function P(x, L,t); that is, the probability that at time t, there is a cell whose length is L and whose center of mass is 
located at x. In Section llVI we use the discrete master equation to derive a partial differential equation for P(x, L, t) 
in a continuous limit which assumes that cell changes its position and length at each Monte Carlo step by a small 
amount. We show that the dependence of P(x, L, t) on L is very close to the Boltzmann distribution. This is used in 
Section^for the derivation of a Fokker-Planck equation for the probability density function p(x, t) of a cell's center of 
mass being at x which is the main result of the paper. In Section IVll it is shown that addition of the time dependence 
of chemical field reduces the Fokker-Planck equation to the Keller-Segel equations. Section lYlII deals with numerical 
verification of the theoretical results of the previous sections and compares the Monte Carlo simulations for our CPM 
and Keller-Segel models. 

II. THE CELLULAR POTTS MODEL 

The Cellular Potts Model (CPM), an extension of the Potts Model from statistical mechanics, is a flexible and 
powerful way to model cellular patterns. Its core mechanism is the competition between the minimization of various 
energy terms in some generalized functional of the cellular configuration, e.g., surface minimization, cell-cell contact 
and chemotactic interactions, and global geometric constraints. It simulates stochastic fluctuations of cell shapes as 
simple thermal fluctuations. 

The CPM is defined on a rectangular lattice £, which is of the form [C^mJ (for 1 dimension), [0,771^] x [0,m y ] (for 
2 dimensions) or [0, m x ] x [0, m y ] x [0, m z ] (for 3 dimensions). (Here [0, m] = {0, 1, ... , m}.) The elements of C are 
called the lattice sites (intervals in ID, pixels in 2D, voxels in 3D). A lattice site is denoted by a index i 6 C. 

Each lattice site has an assigned "spin" <r(i) which can have values s — 0,1,..., Q, where s — corresponds to 
absence of any cell at the given site and the value 1 < s < Q means that the given site is occupied by the sth cell, 
where Q is the total number of cell in the system. Assume that we fix the values of <j(i) at each lattice site, then we 
refer to that set of values as a configuration. The best way to visualize a configuration is to regard the different spins 
as different colors. Each lattice site i has a color cr(i). The cells are the collections of lattice sites that have the same 
spin (color) , so that each lattice site can be occupied by a single cell only. White color correspond to absence of any 
cell at given site: c(i) = 0. In the model considered here we assume that cells cannot divide so that sites with the 
same color are always connected. 

We assume periodic boundary conditions so that pixels at zero position in x, y or z are identical to sites with 
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FIG. 1: (a) Streaming of Dictyostelium discoideum towards the aggregation center. Cells move chemotactically towards the 
aggregation center leading to formation of cell streams and finally mounds. (Reproduced from [2(| with permission), (b) 
Example of a quasi-one-dimensional motion of Dictyostelium discoideum inside a stream (this picture is on much smaller scale 
compared with (a)). Cells are moving parallel to each other in the direction of chemical gradient (from left to right). Chemical 
gradient also causes polarization of cells so that they become elongated in the direction of a gradient. (Reproduced from |lfl ] 
with permission.) (c) Schematic picture of cell motion in a gradient of chemical field (e.g. chemo-attractant cAMP). The 
concentration of the chemical field is shown schematically above the main figure. 



i x = m x + 1, i y = Tfiy + 1 and i z = m z + 1, respectively. 

The temporal dynamics of the system is defined by certain probabilistic transition rules between the configurations, 
giving rise to a Markov chain of configurations, i.e. a sequence of configurations er° , er 1 , er 2 , .... To describe the 
transition rules, we associate to each configuration a an energy E(<r), also referred to as the Hamiltonian. The 
state changes from one configuration to the next are governed by an energy minimization principle with effective 
temperature T. This is implemented by means of the Metropolis algorithm for Monte-Carlo Boltzmann dynamics. 
The algorithm works as follows: 

Given a configuration cr n , we randomly select a lattice site i £ £ such that not all of its nearest lattice neighbors 
have the same spin. We then randomly choose a lattice neighbor i' of i with <j n (\') ^ cr n (i). Let a' be the configuration 
we obtain by "flipping" the spin of i, i.e. we have <r'(j) — cr"(j) for all j ^ i, and cr'(i) = <7™(i'). The new configuration 
<7 n+1 is then either a n or the configuration a'. The probability $(Ai?) that a' is accepted as the next configuration 
cr n+1 depends on the energy difference AE = E(a') — E(a n ). The formula is 

v ; \exp(~(3AE), if AE > w 

Here (3 = 1/T is a positive constant. 

In this paper, we consider a quasi-one-dimensional CPM, which means that cells are assumed to move along x 
direction only and have fixed thickness l y in the y— direction (see Fig. |2J). Let eAa; denote the size of lattice site, 
where < e <C 1, e is the small dimensionless constant and Ax is a dimensional constant of the order of one. Each 
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FIG. 2: Example of a cell in one-dimensional CPM. The cell (shaded domain) occupies lattice sites 2, ... ,6. It has a length of 
5eAx, its center of mass is located at x = 4eAa;, and its end points are xi — 1.5£Aa; and x r = 6.5eAx. 



lattice site is described by its index i = 0, 1, . . . , so that the center of each lattice site is located at x — \eAx with the 
lattice site left border at xi = (i — ^)eAx and the lattice site right border at x r = (i + ^)eAx (see Figure 0) 

In what follows, we will consider the dynamics of a single cell so that the spin a can take two values: if cell is 
absent at a given site and 1 if cell occupies a given site. However, our results remain valid for an ensemble of n cells 
which are well separated from each other, so that the probability that two cells would try to occupy the same volume 
is negligible. This allows us to neglect cell-cell contact interactions. We assume that cells can interact only with the 
substrate (haptotaxis) and the chemical field c(x) (chemotaxis) . The chemical field is assumed to depends only on x 
but not on y. Cells can also produce a chemical which then diffuses. In Section IVD we discuss production of chemicals 
by cells. 

A natural biological realization of this quasi-one-dimensional model is the motion of biolo gica l cells in streams plf. 
E.g. the amoebae Dictyostelium discoideum under starving condition typically forms streams |24j. The biological cells 
inside each stream are moving towards the aggregation center (see Fig. ^i), which results in complicated 2D patterns 
25] . If we zoom to a small scale, we will see that the motion of cells inside each stream is quasi-one-dimensional 
with cells moving parallel to each other in a;— direction ( Fig. \3p). The chemical gradient of the other direction (y 
direction) could be neglected and during cells movement there is no cell-cell interactions, such as cell collisions or cell 
signaling. Fig. schematically shows such a parallel motion of the cells from left to the right under the action of 
the gradient of a chemical field (chemo-attractant). 

For a given configuration a of spins, let N = N(a) denote the number of lattice sites that the cell occupies. The 
length of the cell is equal to L — NsAx. We denote the position of the center of mass of the cell by x and denote the 
position of the left and right ends of the cell by x t and x r , respectively. Then L — x r — x t . (See Figure [3 ) 

We assume that the chemical field c(x) is a slow function of time so its typical time scale is much bigger than the 
time step of a Monte Carlo algorithm. Then the Hamiltonian is given by the formula: 

E = J cm ■ (2L + 2£ y ) + A (£ - L T f + /i c{x) L. (2) 

The first term is a surface energy term which corresponds to the cell-substrate interaction energy (haptotaxis), where 
J cm is an interaction energy between the cell and the medium per unit length. The second term is a length-constraint 
term which penalizes deviations of the cell length L from the target cell length Lt- Here A is a positive constant. 
The choice of A and (i is determined by the typical scale of fiuctiuations of the cellular shape. The third term in © is 
the coupling chemical energy. This term will favor cell motion down or up the chemical gradient for fi > and fi < 0, 
respectively. We assume that the concentration c(x) is a slow function of x on a scale of the typical cell's length L : 

x c /L > 1, (3) 

where x c is a typical scale for variation of c(x) in x. This is consistent with the generally accepted view that cells are 
typically too small to detect chemical gradients without moving. (See e.g. |27| : however recent experimental evidence 
may put this view in question |2^.) Note that the chemical energy could also be defined as /i f Xr c(x)dx. But in the 
limit (PJ), this is equivalent to the form used in the Hamiltonian (|2*|). 



III. DISCRETE EVOLUTION EQUATION FOR PROBABILITY DENSITY FUNCTION 

In this section, we develop an analytical model for the evolution of the stochastic dynamics of a cell in CPM. 
Let P{x,L,t) be a probability density for the cell with the center of mass at x of length L at time t. Spins er(i) 
are defined on the lattice C so that the length of the cell L, which is the difference between positions of right and left 
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ends of cell: L = x r — xi, can take values neAx, n = 1, 2, . . .. The position of the center of mass x = (x r + xf)/2 can 

take values neAx/2, n = 1, 2, That is, the CPM grid is twice the size of the grid of center of mass. In particular, 

if 2^^ is an even number (i.e. x coincides with one of the lattice sites) then the ratio is also an even number. 
Alternatively, if 2j^ is an odd number (i.e. x coincides with a boundary between two neighboring lattice sites) then 
the ratio — k— is an odd number. 

For convenience, we choose a normalization for P(x, L, t) such that the probability for a cell to have its center of 
mass at x and length L at time t is given by (eAx) 2 P(x, L, t). The factor (eAx) 2 results from the product of eAx/2 
(the spacing between lattice sites) and 2eAx (the spacing in L for a fixed x). With this normalization, P(x, L,t) 
becomes a true probability density in the continuous limit e — ► 0. 

We choose the time interval between two Monte Carlo steps to be e 2 At, where At is a fixed constant of dimension 
of time. This implies diffusive time-space scaling, 

e 2 At _ At 
(eAx) 2 = (Ax) 2 

which is independent of the scaling parameter e. We now switch from measuring time in Monte Carlo steps n = 0, 1, . . ., 
to a continuous time variable t = ne 2 At. 

Suppose at time t the cell is at a state (x, L) meaning that it has length L and its center of mass is at x. The 
stochastic discrete system at time t + e 2 At can switch to one of the following four possible states: 

(a) (x + eAx/2, L + eAx) by adding the lattice site x r + eAx to the right end of cell; 

(b) (x + eAx/2, L — eAx) by taking away the site xi from the left end of the cell; 

(c) (x — eAx/2, L + eAx) by adding the lattice site xi + eAx to the left end of cell; 

(d) (x — eAx/2, L — eAx) by taking away the site x r from the right end of the cell. 

Therefore, the most general master equation for evolution of the probability density P(x, L, t) has the form 

P(x, L, t + e 2 At) = [l-T t (x- | Ax, L + eAx; x, L, t) - T r (x + -Ax, L + eAx; x, L, t) 

e e 
—Ti(x + 2^' T ' ^ ~ sAx; x, L, t) — T r (x — t^^ x i L ~ eAx; x, L, tj\ P(x, L, t) 

+Ti(x, L; x + | Ax, L - eAx, t)P(x + -Ax, L - eAx, t) 

+T r (x, L; x - | Ax, L - eAx, t)P(x - | Ax, L - eAx, t) 

+Ti(x, L; x — | Ax, L + eAx, t)P(x - -As, L + eAx, t) 

+T r (x, L; x + | Ax, L + eAx, t)P(x + ~Ax, L + eAx, t), (4) 

where Ti(x, L; x', L') and T r (x, L; x' , L') correspond to transitional probabilities for a cell of length V and center of 
mass at x' to change into a cell of length L and center of mass at x'. Subscripts "1" and "r" corresponds to transition 
due to addition/removal of a pixel from the left/right side of a cell respectively. These transition probabilities are 
given by 

Ti(x,L;x',L')=T r (x,L;x',L') - ^(e(x,L) - E(x' ,L')^j, (5) 

where E(x,L) is the Hamiltonian J5J and $>(AE) is given by Eq. 0). Factor 1/4 in J^J accounts transitions to 4 
possible states (a)-(d). For computational purposes it is convenient to rewrite Q in an equivalent form 

$(AE) = l-|l-exp[-/3A£;]}e(AJ5). (6) 

Here O(x) is a Heaviside step function: 0(x) = 1 for x > and O(x) = for x < 0. 
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IV. CONTINUOUS EVOLUTION EQUATION FOR PROBABILITY DENSITY FUNCTION OF CPM 

Below we assume e to be small, e <C 1, so that the change of the cell size and position is small at each Monte-Carlo 
step. Now we carry out a Taylor series expansion in e of the terms in Eq. (0J. One has to take special care of 
0(A_E) terms in the expansion because the Heavyside step function is not analytic. To avoid this difficulty we do 
not expand the function itself but only its argument instead. There is an important simplification which comes from 
the fact that Q(AE) + Q(-AE) = 1 so that in Eq. gj| we obtain that T l>r (x,L; x',L',t) + T iir (x' ) L'; x,L,t) 

1 1 . exp ( — (3\E(x, L) — E(x' , L')\) . This yields mutual cancellation of nonanalytical terms up to order 0(e 2 ). Then, 
equating coefficients in the Taylor expansion in Eq. in order 0(e 2 ) results in the Fokker-Planck equation 

d t P(x, L, t) = D{d 2 x + 4d 2 L )P + %D(3Xd L {LP) + D(3L^d x [c{x)P)} , 

L = ~[j cm + \(L-L T ) + ~»c(x)], = ^^. (7) 

Now, under certain conditions to be described in the end of this section, the terms DAdj^P + &D(3XdL{LP) dominate 
the other terms on the right hand side of Eq. JJJ. This means that at the leading order, one can neglect terms with 
x— derivatives. Under this assumption, the probability density function P(x, L, t) approaches a Boltzmann distribution 
for cell length exponentially in time at the rate of 8Df3X: 

P(x, L, t) = P B0 Uz(x, L)p(x, t), (8) 

where p(x, f) is a probability density function of finding cell's center of mass at x. PboUz{x,L) is the Boltzmann 
distribution for the cell length given by 

PBoitz{x,L) = i exp(-PAE lengt h), (9) 

^Elength = E(L) — E mm = XL 2 , (10) 

where E m i n is a minimum of energy E(L) as a function of L for a given x, 

E m in — E(L m i n ), L rn i n — Lt — — — , (11) 

and Z is a partition function 

Z(x) = 2eAx x ^ exp(-PAEien g th), 



L=(l+a)sAx, (3+a)eAx, (5+a)eAx, 

X nt x 
— — = n, a = (J tor — — 

sAx eAx 



a = 1 for — — = n, a = for = n + 1/2, n G N. (12) 



Here we use the fact that due to discrete nature of our model, the position of the center of mass, x, could be located 
at one of the lattice sites x = me Ax (m being an integer number) if the length of the cell L is an even number of 
units eAx or x could be located at the boundary between two neighboring lattice sites in case of L being equal to 
an odd number of units of eAx. The factor (eAx) 2 in the definition of the partition function H12(l is chosen in such 
a way as to yield J P(x, L,t)dLdx = 1 in the continuous limit. We can also normalize J P(x, L,t)dLdx — N to the 
total number of cells in the system N. 

In the continuous limit, e — > 0, the sum in Eq. i|12|) is transformed into the integral 

exp(-0AE lmgth )dL=^=, x^O. (13) 
-oo VPA 

Here we have extended the limits of integration from (0, +oo) to (— oo, +00). Of course physically, the length of the cell 
L is always positive. A typical fluctuation of the cell size 8L = L — L m i n about L m i n is determined by the Boltzmann 
distribution JSJ as (3X8 L 2 ~ 1. In what follows we make a biologically motivated assumption about fluctuations of 
the cell size being much smaller than L: \8L\ <C E m in which results in the condition 

pL 2 min X » 1. (14) 

This justifies the use of the integration limits (—00, +00) in Eq. (|13fl instead of (0, +00) because under this condition 
exp(— (3AEi eng th) peaks around L m i n and replacement of integration limits results in an exponentially small correction. 
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Let us now specify the conditions for the applicability of the Boltzmann distribution approximation JSJ). For this, 
consider Eq. (JJJ. We have (3\5L 2 ~ 1. We now assume in addition the relation 

/3^A>1, (15) 

where xq is a typical scale of P with respect of x. Note that under the assumption that L m i n <C xq, i.e. that 
the typical length of a cell is much smaller than xq, the condition (|15() follows from (|14|l . It follows from 115(1 that 
\d%P\ *C |49fP|, and consequently we may neglect the first term with x— derivative , d xx P, on the right hand side of 
Eq. 0. 

The second condition for the applicability of the Boltzman distribution approximation (jSJ) is the assumption that 
the last term with x— derivative in Eq. Q is small, \{3Lfj,d x [c'(x)P] | <C |4d£P|. This is true if 

\L m inVCo\(l + — ) < Ax 2 ,, (16) 

where Co is a typical amplitude of c(x) and x c is a typical scale of variation of c(x) with respect to x. Lastly, recall 
that we derive the continuous Eq. (JJJ from the master equation Q under the condition of the step in x being small 

e«l. (17) 

Notice that diffusion coefficient D in Eq. Q does not depend on j3. Instead j3 determines a rate of convergence 
TjT 1 = 8D/3X of P(x,L,t) to the Boltzmann distribution |jSJ. 

We have solved both the master equation (J3J and its continuous limit J7J) numerically with initial conditions 
P(x,L,0) different from the Boltzmann distribution 10}. Simulations described in Section IVlII demonstrate that for 
each x, the solution P(x,L 1 t) indeed converges in time to the Boltzmann distribution at an exponential rate of 
- 8D(3\. 



V. FOKKER-PLANCK EQUATION FOR PROBABILITY DENSITY FUNCTION p(x,t) 



We now turn to calculating the probability density function p(x, t) of a center of cell's mass being at x. It is given 
by the sum over all possible lengths of a cell 



p(x, t) = 2s Ax 

L = (l+a)eAi, (3+a)eAi, (5+ajeAi, 

nr 

a = 1 for 



sAx 



/-too 
P(x,L,t)dL, e->0, 
-CO 

X 

n, a = for = n + 1/2, n G N, 



eAx 



(18) 



which reduces to Eq. (jHJ in the Boltzmann distribution approximation limit. 

To derive closed equation for p(x, t) we substitute ansatz JSJ into and integrate both right hand and left hand 
sides of Eq. Q with respect to L to obtain 



X(x) = —0fi J cm - \L T + ^fic(x) 



d t p = Dd 2 x p - d x [x{x)p d x c(x)] , 

(Ax) 2 



D 



lAt 



(19) 



This continuous equation is the main result of this paper. The conditions for the applicability of Eq. (|19|l are given 
by Eqs. flU, C3> 03 and $TQ. 



VI. REDUCTION TO KELLER-SEGEL MODEL 



In this section we add time dependence to the chemical field c (concentration of chemoattractant or chcmorepellant) 
by including a diffusion equation with the source term ap which determines the secretion of chemical by a cell 

d t c = D c d 2 x c- 1 c + ap 1 (20) 

where D c is a diffusion coefficient of the chemical field, 7 is the decay rate of the chemical field and a is a production 
rate of the chemical field. 
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The system of equations <|19f) and l|20|) is applicable under the assumption that the typical time scale r c of diffusion of 



c(x, t), given by r c = ^ SL , is large in comparison with convergence time r r — 1/(813 /3A) of P(x, L, t) to the Boltzmann 
distribution JHJl, where x c is a typical spacial width of the distribution of c(x, i). Namely, this condition has the form 

T c /r r = 8D/3Xt c > 1, (21) 

Eqs. I|19|) an< i l|20() form a closed set of equations which is equivalent to the classical Keller-Segel model p] of 
chemotaxis. If the parameters satisfy condition 

\J cm -XL T \>^\»\c(x), (22) 

than Eq. (|19Jl reduces to the following commonly used form of the Keller-Segel model 0, El : 

d t p = Dd 2 x p - xad x [p d x c] , 
Xo = DXPfi [,J cm - XLt\ , D = . (23) 

The probability density function p(x, t) corresponds to the microscopic density in the Keller-Segel model. Notice that 
both in the Keller-Segel model and CPM considered in this paper, there is no direct interaction between cells except 
through production and reaction to a chemoattractant. In other words, cells are treated in a way similar to a dilute 
gas with long range nonlocal interactions due to reaction to a chemical field. 



VII. COMPARISON OF NUMERICAL SIMULATIONS 



In this section, we describe numerical tests comparing Monte Carlo simulations of the CPM and simulations of both 
discrete and continuous models for the probability density functions P(x,L,t) and p{x,t), as given by Eqs. (JTJ) 
and dl . 



A. Monte Carlo simulations 



The computation of the frequency distribution of the cell center of mass and length for the CPM has been carried 
out as follows: 

1. We run a large number N of CPM simulations with one cell with the same initial conditions. 

2. We fix a time interval St = e 2 At, i.e. we fix the time interval between successive Monte Carlo steps. For each 
simulation we record the locations of the center of mass and and lengths of the cell at the times t = St, 2St, 3St, . . . 

3. After the N runs, the recorded data give a frequency distribution M(x,L,t) for the location of the center of 
mass of the cell and length of the cell. 

The frequency distribution M(x,L,t) determines the approximation P cpm (x, L,t) — M(x, L,t)/(N(eAx) 2 ) of the 
probability density function P(x,L,t) for the center of mass of a cell of length L being at x at time t. Therefore, 
we compare P cpm (x, L,t) with P(x,L,t) which is a solution of either the master equation Q or the Fokker-Planck 
equation J7J). To approximate the probability density function of center of mass p(x, t) we sum up over all values of 
L on the grid in a way used in Eq. (|18|) 

Pcpm(%jt^ — 2s Ax ^ ^ Pcpm (■£ , t) , 

L=(1+q)eAi, (3+a)eAx, (5+a)eAi,... 
X X 

a = 1 for — — = n, a — for — — = n +1/2, n e N. (24) 
eAx eAx 

In what follows, we compare p C pm(£, t) for s 1 with p(x, t), a solution of the continuous Eq. I|19f) . corresponding to 
the following choice of parameters 



A = 4, L T = 5, J cm = 2, (3 = 15, /x = 0.1, Ax = 1, At = 1. 



(25) 
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The size of the CPM lattice is chosen to be L cpm = 100; and the model is typically run from to = to t en d = 200. The 
number of the CPM lattice sites and the number of Monte Carlo steps are chosen to be L °££ and Jf^r respectively. 
We use a range of values of e between 0.2 and 0.001. 

The initial conditions for each CPM run are chosen as follows. A random pixel in the interval [40, 60] is selected as 
a center of mass of a cell, and then the length L for the cell is chosen with probability Z^ 1 exp(E(L) — E(Lt))- Here 
the normalization constant Zi is chosen to have the total probability 1. In most simulations, we use the following 
distribution for the chemical field c{x): 



B. Monte Carlo simulations versus numerical solutions of the discrete master equation and the 

Fokker-Planck equations 

We first compare Monte Carlo simulations with the numerics for the master equation Q and the Fokker-Planck 
equation (7J . Simulations of the Fokker-Planck equation 0) have been performed by using a finite-differences scheme. 
Figure |3 shows the probability density functions for all three types of simulations. The difference between the master 
equation an d the Fokker-Planck equation Q simulations is negligibly small for e = 0.01 (Fig. [3^) but can be 
clearly seen for s = 0.1 (Fig. EJs). We conclude that for N — > oo, the Monte Carlo simulations converge to the solution 
of the master equation (@J for any e. The rate of convergence is about N^ 1 / 2 . For small e — > 0, the solution of the 
Fokker-Planck equation JJJ also converges to the solution of the master equation. 



C. Convergence of the probability density function P(x,L,t) to the Boltzmann distribution 

To demonstrate quick convergence of P(x,L,t) to the Boltzmann distribution (JHJ) (as discussed in Section Hvj) 
we solve numerically both the master equation @ and its continuous limit J7J with initial conditions P(x, L, 0) 
being different from the Boltzmann distribution jljl. Namely, we choose initial value P(x, L,0) to be the Bolzmann 
distribution with different temperature 0i n i = 1.5 so that Figure^Jshows convergence of initial state with temperature 
1/ ' f3i n i to the quasi-equilibrium state with temperature 1//3 = 1/15 used in the Monte Carlo algorithm. Linear-log plot 
in the Figure ^indicates that convergence is indeed exponential in time with high convergence rate t~ x (t^T 1 = 98.55 
for parameters of Fig. • By high convergence rate we mean that the typical convergence time r r is small compare 
with e.g. the diffusion time Xq/D in x (see Eq. (JJJ). Because of the x-dependence of the chemical field, the convergence 
rate r r is also x-dependent and a closed analytic expression for it is difficult to obtain from Eq. Q for general c(x). 
However even a simple estimate r" 1 = 8D/3X of the rate of convergence gives 60 for parameters of Figure 0] which is 
qualitatively close to numerical value 98.55. Here 98.55 is obtained from the linear fit presented in Figure^] 

Also, we observe that if we increase temperature T in Monte Carlo simulations, so that condition 1|14[> is not true 
any more, then it results in a significant departure from the Boltzmann distribution (jSJ which confirms the theoretical 
results of Section ITVI 



D. P(x,L,t) vs. p(x,t) simulations 

The ansatz (jSJ can be used for fast simulations of solutions of the discrete master equation. Summing up over all 
values of L in the master Eq. JIJ and taking into account (JSJ result in a discrete equation for the probability density 
function p(x, t) 

p(x, t + e 2 At) = [l - T(x + ~Ax; x, t) - T(x - |Ax; x, t)]p(x, t) 
+T(x\ x - |Ax, t)p(x - ^Ax, t) + T{x; x + -Ax, t)p(x + |Ax, t), (27) 

where T(x; x' , t) is a transition probability of a change of position of a center mass from x' to x at time t. Expressions 
for T(x] x' , t) arc described in the Appendix. They are calculated only once at the beginning of a simulation which 
makes the numerics for discrete Eq. (I27|l very efficient. 

We run simulations for the discrete equation l|27ll and the continuous Eq. (|19l) and compared them with the solutions 
of the discrete (@J and continuous (JZJ) equations, respectively. We find, taking into accout Eq. (jHJ, that indeed the 
differences between these solutions are very small for the typical values of parameters. 
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FIG. 3: Probability densities for Monte Carlo simulations p cp m(x,t) 
(dotted line), p(x,t) for the Master Eq. (solid line) and the Fokker-Planck equation (0 (dashed line) versus x for t = t en d- 
(a) e = 0.01; (b) e = 0.1. The difference between position of solid curve and a dashed curve is negligibly small in (a). Number 
of Monte Carlo simulations is N = 2 x 10 . We used c(x) as given by l|26|l . 

We conclude that the Monte Carlo simulations of CPM are equivalent in the limit of large N to the the simulations 
of the discrete Eq. it27|) for any e. 

E. Comparison of the continuous model with the CPM 

Below we denote as p cpm both Monte Carlo simulations and numerical solutions of Eq. (|27|l and as p CO nt{x,t) 
solutions of l|19fl . 
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FIG. 4: Exponential convergence of the Full Width Half Maximum (FWHM) of P(x, L, t) in L as a function of time for x = 50. 
The vertical axis corresponds to the normalized difference [VJ^(t) — W@]/W@, where W(t) is the FWHM at time t and Wp is 
the FWHM for the Boltzmann distribution (29 . Solid squares correspond to the numerical solution of both Eqs. IgJ and Q. 
The solid line is the best linear fit which gives exponential convergence e _98 55t . The same parameters as in Figure[3]are used 
here with e = 0.01. 



Figure [SI shows a series of simulations of the CPM (dotted line) and numerical solutions of the continuous Eq. 1191) 
(solid line) for different values of e. This Figure demonstrates that in the limit e —> 0, the solution of the continuous 
Eq. Ijl 9(1 appears to converge to the cell probability density function of the CPM. 

Figure El shows the normalized difference between solutions of ((19(1 and the CPM. The normalized difference ap- 
proaches as e decreases. We also run a series of tests for different forms of the chemical field c(x) and demonstrate 
that solutions of the CPM and continuous Eq. ((19(1 are close for small values of e. Figure [7] shows a typical result 
of numerical simulations for a "double well" chemical concentration c(x) = cos(47ra;/100). We conclude that the nu- 
merical simulations show excellent agreement between the CPM and the continuous Eq. 1 (19(1 provided that the Potts 
parameters satisfy conditions 1(14(1 . 1 (15(1 . 1(16(1 . ((17(1 and e — > 0, which correspond to the continuous limit of the CPM. 

VIII. CONCLUSIONS 

In this paper we combine microscopic and macroscopic levels of description of one dimensional cellular dynamics. 
The microscopic level is represented by a one dimensional CPM with chemotaxis and without cell-cell adhesion 
term. We study a continuous macroscopic limit of our CPM as the size of Monte-Carlo step is made small under 
the assumption that changes in the cell's position and length are also small. In this limit, we derive the Fokker- 
Planck equation 1(19(1 for the probability density function p(x, t) of cells and then further reduce it to the well-known 
macroscopic continuous Keller-Segel model 1(20(1 and 12: il) for the chemotactic aggregation of cells. All coefficients of 
the Keller-Segel model are derived from parameters of the CPM. 

We use numerical simulations to test hierarchy of models and assumptions which we used to derive continuous 
equation ((19(1 . In particular, we compare Monte Carlo simulations with simulations of both the discrete master 
equation Q and the Fokker-Planck equation Q for P(x,L,t). We find that, as expected from our theoretical 
analysis, all models agree for small e. Also Monte Carlo simulations agree with the solutions of the discrete master 
equation (3J for arbitrary s. We verify numerically that the probability density function P(x,L 7 t) quickly converges 
to the Boltzmann distribution (5J . And finally, we find that numerical simulations show excellent agreement between 
Monte Carlo simulations of CPM and the continuous macroscopic model ((1911 . 

We are currently working on extending our results to a 2D case for modeling chondrogenic patterning in the presence 
of chemotaxis and fibronactin production (2t| . 
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FIG. 5: Plots of Pcpm (dotted line) and p CO nt(x,t) (solid line) as functions of x for a series of decreasing values of e at time 
t — 200. All other parameters are the same as in Figure 
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FIG. 6: Normalized difference between solution of CPM and continuous Eq. I|19|l for the same parameters as in Figure as a 
function of e. Normalized difference is given by 1 — f p cpm (x,t)p con t(x,t)dx/ J p cont (x,t) 2 dx for t = t an d- 
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FIG. 7: Typical results of CPM simulations. The same parameters as in Figure [3] are used except that c(x) — cos(47ra;/100), 
e — 0.01. The same notation for solid, dashed and dotted curves as in Figure is used here. The difference between position 
of solid curve and a dashed curve is again negligibly small. 



X. APPENDIX 



The explicit expressions for the transitional probabilities T{x\ x' , t) used in Eq. (|27|l can be obtained by summing 
over all lengths (or, in other words, over even multiples of eAx (if 2x/(eAx) is an even number), and over odd 
multiples of sAx (if 2x/(eAx) is an odd number). A change in the position of the center of mass from x to x ± |Aa; 
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can be made by adding/removing lattice sites from the left/right end of a cell which results in 

T(x; x - -Ax, t) = AZ(x _ E _ Ax) E 

| exp [ - fiAEiength (x - ^Ax,L- eAx)]$(e(x,L) - E(x - ~Ax,L- sAx) 
+ exp [ - pAEi en gth {x~^Ax,L + eAx)] $ (e(x, L) - E(x - | Ax, L + eAxfj } 



T(*;x+-A M ) = 4 E 

^ ' L=(1+q)eAi, (3+a)eAi, (5+a)EAi,. 



| exp [ - (3AE length (x+^Ax,L- eAx)] $ (e(x, L) - E(x + ^Ax,L- eAx) 
+ exp [ - /3AEi ength (x + | Ax, L + eAx)] $ (i?(x, L) - B(a; + | Ax, L + eAx)) | 
T(x +- Ax; x, <)-—-- E 

v ; t=(l+a)EAi, (3+a)eAi, (5+o)eAi,.. 

| exp [ - f3AEi ength (x, L)] $ (e(x + |Ax, L - eAx) - £(x, L) 
+ exp [ - 0AE length (x, L)] $ (£(x + | Ax, L + eAx) - £(x, L)j j 

T(x--Ax;x,t) = i ^y E 

v ; L=(1+(i)eAi, (3+q)eA3:, (5+o)eAi,.. 

{ exp [ - {3AE length (x, L)] $ {e{x - |Ax, L - eAx) - £(x, L) 
+ exp[- 0AE length (x,L)]^{E(x- ^Ax, L + eAx) - E(x, L)j j, 

a = 1 for -4— = n, a = for = n + 1/2, neN. (28) 
eAx eAx 

Here the partition function Z(x) is given by l|12fl . Z(x) is x-dependent in the discrete case considered in this Appendix. 
This x-dependence is eliminated after going from a discrete summation in (|12fl to an integral (as in Eq. I|13|) h We 
evaluate the transitional probabilities T(x; x ± § Ax, t) and T(x ± § Ax, t) numerically using (|28(1 for each value of x 
once at the beginning of each simulation and then calculate the discrete evolution of Eq. Ij27(l . 

Notice that in the limit of small e — > 0, the continuous equation l|19l) can be derived directly from l|12l) . <|27[1 and 
(|28|l . However, this derivation is more tedious compared with the two-step derivation in Sections IIVI and Ivl where 
continuous equation JJJ) is first derived and then integrated (0) over L which results in Eq. (|19|l . 
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